// Function to setup simulation 
void initSimulation() {

	ifstream simInput; 
	
	simInput.open("rev_micelle_mc.inp"); 

	if (!simInput) { 
		cerr << "File Open Error" << endl; 
		exit(1); 
	} 

	else { 
		cout << "Reading simulation parameters ....." << endl << endl; 
	} 

	char comm1[50];
	simInput >> comm1; 

	simInput >> deviceID >> nMolecules >> nSteps; 
	simInput >> comm1; 
	
	simInput >> tempK >> atomPressure >> simulationType >> nptSelect; 
	simInput >> comm1; 

	cout << nptSelect; 

	simInput >> moveSize >> dRm >> deltaLmax; 
	simInput >> comm1; 

	simInput >> SIGMA_OO >> SIGMA_CC >> SIGMA_OC; 
	simInput >> comm1; 
	
	simInput >> EPS_OO >> EPS_CC >> EPS_OC; 
	simInput >> comm1; 
	
	cout << "Using GPU # " << deviceID << endl << endl;
	cout << "SIMULATION PARAMETERS: " << endl << "-----------------------" << endl; 
	cout << "Running " << simulationType << " for a system of " << nMolecules << " CO2 molecules for " << nSteps << " steps" << endl << endl;
	cout << "Temp (K): " << tempK << "\t | Pressure (Pa): " << atomPressure << endl << endl;
	cout << "Move Size: " << moveSize << "\t | Rotational Move Size: " << dRm << "\t | Max. Box Change: " << deltaLmax << endl << endl; 
	cout << "CO2 MODEL PARAMETERS: " << endl << "----------------------" << endl; 
	cout << SIGMA_OO << "(sigma_oo)\t" << SIGMA_CC << "(sigma_cc)\t" <<SIGMA_OC << "(sigma_oc)\t" << endl;
	cout << EPS_OO << "(epsilon_oo)\t" << EPS_CC << "(epsilon_cc)\t" << EPS_OC << "(epsilon_oc)\t" << endl << endl;

} 

// Function to read positions 
void readConfiguration() { 

	int total; 
	char type; 

	// Equilibration: Read from the initial configuration file //
	if (simulationType == "EQ") { 
		
		ifstream eqInput; 
		eqInput.open("active_configuration.xyz"); 

		if (!eqInput) {
			cerr << "File Open Error" << endl;
			exit(1);
		} else {
			cout << "Reading reference configuration ..... \n"; 
		} 
	
		eqInput >> box.x >> box.y >> box.z; 
		eqInput >> total; 
          
		for (int j=0; j<nMolecules; j++) {
			eqInput >> type >> hAtom[j*3+2].x >> hAtom[j*3+2].y >> hAtom[j*3+2].z; 
			eqInput >> type >> hAtom[j*3+1].x >> hAtom[j*3+1].y >> hAtom[j*3+1].z;
			eqInput >> type >> hAtom[j*3+0].x >> hAtom[j*3+0].y >> hAtom[j*3+0].z;         
		}
                
	}       

	// Production: Read from the equilibrated configuration file //
	if (simulationType == "PROD") { 

		ifstream prodInput; 
		
		prodInput.open("eqconf.out"); 
		prodInput >> box.x >> box.y >> box.z; 
		prodInput >> moveSize; 
		prodInput >> dRm; 
		prodInput >> deltaLmax; 

		for(int j=0; j<nMolecules; j++) {

			prodInput >> type >> hAtom[j*3+2].x >> hAtom[j*3+2].y >> hAtom[j*3+2].z;
			prodInput >> type >> hAtom[j*3+1].x >> hAtom[j*3+1].y >> hAtom[j*3+1].z;
			prodInput >> type >> hAtom[j*3+0].x >> hAtom[j*3+0].y >> hAtom[j*3+0].z;   
		}
	}
} 

// Function to read positions 
void writeConfiguration() { 

	if (simulationType == "EQ") { 

		ofstream eqLog;

		eqLog.open("eqconf.out"); 
		eqLog << box.x << "\t" << box.y << "\t" << box.z << endl;
		eqLog << moveSize << endl << dRm << endl << deltaLmax << endl;

		for(int i = 0; i<nMolecules; i++) {
	
			eqLog.precision(10); 
			eqLog << "C\t" << hAtom[i*3+2].x << "\t\t" << hAtom[i*3+2].y << "\t\t" << hAtom[i*3+2].z << endl;
			eqLog << "O\t" << hAtom[i*3+1].x << "\t\t" << hAtom[i*3+1].y << "\t\t" << hAtom[i*3+1].z << endl;
			eqLog << "O\t" << hAtom[i*3+0].x << "\t\t" << hAtom[i*3+0].y << "\t\t" << hAtom[i*3+0].z << endl;

		}
	} 

	if (simulationType == "PROD") { 
		ofstream prodLog; 
		prodLog.open("prodtraj.xyz", ios::app); 	
		prodLog << nMolecules*3 << endl << "CO2" << endl; // << box.x << "\t" << box.y << "\t" << box.z << endl; 
		for(int i = 0; i<nMolecules; i++) {
			prodLog << "C\t" << hAtom[i*3+2].x << "\t\t" << hAtom[i*3+2].y << "\t\t" << hAtom[i*3+2].z << endl;
			prodLog << "O\t" << hAtom[i*3+1].x << "\t\t" << hAtom[i*3+1].y << "\t\t" << hAtom[i*3+1].z << endl;
			prodLog << "O\t" << hAtom[i*3+0].x << "\t\t" << hAtom[i*3+0].y << "\t\t" << hAtom[i*3+0].z << endl; 
		}	
	}
}
